Gauss collocation methods for efficient structure preserving integration of 

post-Newtonian equations of motion 
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In this work, we present the hitherto most efficient and accurate method for the numerical inte- 
gration of post-Newtonian equations of motion. We first transform the Poisson system as given by 
the post-Newtonian approximation to canonically symplectic form. Then we apply Gauss Runge- 
Kutta schemes to numerically integrate the resulting equations. This yields a convenient method for 
the structure preserving long-time integration of post-Newtonian equations of motion. In extensive 
numerical experiments, this approach turns out to be faster and more accurate i) than previously 
proposed structure preserving splitting schemes and ii) than standard explicit Runge-Kutta meth- 
ods. 
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INTRODUCTION 



When Einstein gave birth to general relativity with 
the presentation of his field equations in 1915, new phe- 
nomena such as black holes and gravitational waves were 
soon predicted as consequences of this theory. In the last 
couple of years, gravitational waves have attracted ever 
more attention. With the aim to finally receive signals 
of such waves, much experimental effort has been put 
upon mounting land-based detectors. Virgo in France 
and Italy, GEO 600 in Germany and the UK, and LIGO 
in the USA are only to name a few. They are soon to be 
joined by the space-based eLISA. In order to track any 
signal of gravitational waves, templates are required that 
give a hint on which needle to look for in the haystack of 
data delivered by all the working detectors. Such tem- 
plates, in turn, can only be obtained by singling out the 
most promising sources of gravitational waves and calcu- 
lating their motion in phase space. The main source of 
waves have been identified to be binary systems consist- 
ing of inspiraling compact objects, see, e.g., [1]. Their 
mass proportions can be anything between equal masses 
and extreme ratios. Binaries with very unequal masses 
are called Extreme Mass Ratio Inspirals (EMRIs). One 
common example of an EMRI is a neutron star that or- 
bits a super massive black hole (SMBH). EMRIs allow 
for a simple description as a free particle (the lighter 
one) moving in a curved spacetime given by the metric 
corresponding to the mass of the heavier particle. Many 
possible shapes of the background metric have been pro- 
posed in this field, e.g., [2, 3]. 

Binaries with not so extreme a mass ratio arc suit- 
ably described by the post-Newtonian formalism. This 
approach was possible after Arnowitt, Deser and Misner 
discovered that Einstein's theory can be formulated as 
a Hamiltonian System, [4]. The idea is then to expand 
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the elements of the metric tensor and the equations of 
motion of the matter in powers of the small parameter 
■i^, see, e.g., [5]. This gives the Hamiltonian as a power 
series in the small parameter, the first term of the series 
being the Hamiltonian for Newton's law of gravitation. 
The determination of the individual terms in this expan- 
sion is subject to current research in theoretical physics 
and contributions up to 3PN order have been given in [6] . 
The post-Newtonian approach has even been extended to 
a binary which is perturbed by a much lighter third body, 
e.g., [7]. 

A property of relativistic test-particles which is not 
known from classical mechanics is their spin. After the 
foundation for the treatment of this spin had been laid 
down in the 1950s, e.g., [8], the post-Newtonian formal- 
ism could be expanded to include the corresponding con- 
tributions. These comprise spin-orbit as well as spin-spin 
interactions, see, e.g., [9, 10]. With this extension, the 
Hamiltonian system becomes a so called Poisson system. 

One important property of the post-Newtonian is that 
they are generally non-integrable. As a consequence, the 
motion described by them can exhibit chaotic traits. If 
the motion of a particular binary is chaotic, the grav- 
itational waves emitted during its inspiral will be un- 
predictable, thus leaving the researchers at the various 
wave detectors without any useful template. Hence, the 
investigation for chaos of a given binary system is an 
important task. Consequently, many works have been 
published concerning this topic both in the geodesic and 
the post-Newtonian field, e.g., [11-15]. The analysis of 
chaos requires reliable indicators and, above all, numer- 
ical simulations over very long time spans. Numerical 
long-term analysis, in turn, relies on there being efficient 
and highly accurate integration schemes which behave 
well even during long-time simulations. To this aim one 
can make use of the post-Newtonian equations' special 
structure. 

Over the last few decades, the numerical analysis com- 
munity came up with tools for the long-term integra- 
tion of equations of motion. In the course of this, struc- 
ture preserving algorithms such as symplectic schemes 
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for Hamiltonian systems (e.g., [16-18]) or symmetric in- 
tegrators for time- reversible systems (e.g., [19, 20]) have 
been proposed. Regarding long-time behaviour and con- 
servation properties, these schemes are superior to ordi- 
nary numerical integrators such as exphcit Runge-Kutta 
schemes in many apphcations of classical mechanics and 
astronomy. Whereas for standard integration schemes 
the overall error is normally proportional to the square 
of the length of the integration interval U, it only in- 
creases linearly with ti for structure preserving integra- 
tors. And whereas there is a drift in constants of motions 
for standard methods, these constants are conserved up 
to a small error over extremely long times for symplec- 
tic algorithms. These algorithms have been successfully 
applied even in quantum mechanics. A comprehensive 
presentation of such methods is given by [21]. 

In the last years, attempts have been made to con- 
struct structure preserving integrators for relativistic sys- 
tems of compact binaries. Most recently, such a scheme 
has been proposed for the geodesic approximation, see 
[22]. In the realm of post-Newtonian equations, two dif- 
ferent approaches have been considered so far. First, 
a non-canonically symplectic integrator has been con- 
structed which preserves the system's Poisson structure, 
see [25]. Then, a transformation to canonical form has 
been proposed, sec [23], after which symplectic meth- 
ods have been applied, [24] . All the previous approaches 
have in common that they are based on a splitting of the 
Hamiltonian into a Newtonian part and other relativis- 
tic contributions. In this work, we will first argue and 
then demonstrate via numerous experiments that a more 
efficient and accurate method for the solution of post- 
Newtonian equations consists of a transformation to sym- 
plectic form followed by the application of Gauss Runge- 
Kutta schemes. This will drastically reduce the numeri- 
cal effort when simulating post-Newtonian systems. 

Our work is organized as follows: We first explain 
our notation in Section 2. Afterwards we will discuss 
the post-Newtonian equations of motion and their nu- 
merical properties in Section 3. In Section 4 we briefly 
summarize the main aspects of the Poisson integrator of 
[25]. Section 5 deals with the transformation to canon- 
ical form. The subsequent Section 6 presents common 
splitting methods. Then, in Section 7 we present Gauss 
Runge-Kutta schemes and argue why they are a good 
choice in post-Newtonian simulations. Finally, we sub- 
ject the individual methods to extensive tests and com- 
pare them to standard explicit schemes in Section 8 be- 
fore we summarize our main results in Section 9. 



notation we will combine the relevant variables of the 
phase space into one variable y = (p,x. Si, 82)"^. With 
this abbreviation we can write the equations of motion 

as 



dy 

dt 

y(0) =yo, 



(1) 
(2) 



with an appropriate function /. An exact solution of 
a differential equation of the form (1) which starts at a 
given point yo and propagates the system over a time t 
will be denoted by 



y(*) = 'Pt(yo), 



(3) 



whereas a numerical approximated fiow over a time step 
h is written as $/i(yo). Consequently, for a given point 
of the phase space y„, the next point on the numerical 
trajectory is calculated as 



(4) 



Finally, J is a unit matrix of appropriate dimension and 
J is the symplecticity matrix 



/ 

-/ 



(5) 



POST-NEWTONIAN EQUATIONS OF 
MOTION 



We consider orbital contributions up to order 3PN as 
given in [6] and the leading term of the spin-orbit (SO) 
and the spin-spin (SS) contribution of [10], respectively. 
With all these terms, our Hamiltonian reads 



if (p,x, 81,82 



=iiorb(P,x) -f-i/so(P,X, 

-^iJss(p,x,Si,82). 



81,82) 



(6) 



2. NOTATION 

In this work we use canonically conjugate position and 
momentum variables and restrict ourselves to the center- 
of-mass frame so that we have the relevant variables (with 
a = 1,2 denoting the individual compact objects) x = 
X2— xi, p = Pi = — P2, 81 and 82. For the sake of shorter 



As the motion docs not depend on the absolute value 
of the masses but on their ratio, one can without loss 
of generality assume the total mass m := mi -t- 7712 to be 
equal to 1. Using the reduced mass /i := ™]^"'^ , q := ||x||, 
1/ := ^, and the unit vector n := ^ and choosing units 
such that G = c = 1. the relevant terms of the orbital 
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Hamiltonian i^orb are 



2pt q 

i^lPN,orb(p,x) = ^(^^^ 1)(P^ 



(7) 



2PN,or: 



b(p,x) 



1 



[(3 + i^)p2 + i.(n . p)2] 



2g2 



(8) 



1 



+ ^[(5-20z.-3^.2)(p2)2- 
2-i.2(n.p)2p2+3,.2(n.p)4] 

+ -i^[3t.(n.p)2 + (5 + 8i/)p2] 
2/i(72 

(l + 3i')/i 



J^3PN,orb(P, x) 



4(73 
1 



(9) 



2\4 



128Ai 
1 



y(-5 + 35i^ - 70j/2 + 35j/3)(p2) 
[(-7 + A2v - 53J/2 - 5i/3)(p2)a 



16^^(7 

+ (2- 3^.^2(11. p)2(p2)2 
+3(1 - ■ p)V - 5i^3(n . p)6j 

1 

+ (17 + 30j/)j/(n-p)2p2 



2\2 



(-27+ 136zy+ 109z/2)(p2) 



-(5 + 43i/)i^(n-p) 



1 



25 



7r2 335 \ 23 . 



64 48 / 8 



85 37r2 7u 



/>4 



16 64 4 , 
"1 /109 2l7r2 
12 32~ 



(10) 



The leading order spin-orbit coupling can be expressed 
by means of the orbital angular momentum L = x x p 
and the effective spin 



Scff = 1 



3m2 
4m 1 



Si 



1 + ^')S2 

47715 



(11) 



The spin-spin interaction is the sum of the three following 
terms: 



Hs,s, (P, X, Si, S2) = ^ [3(Si • n)(S2 • n) - Si • S2] 



i?SiSi(p,x,Si) 

-H"S2S2(P>X, S2) 



7772 



2777ig3 

7771 
277729 



[3(Si-n)2-Si-Si] 
3 [3(S2-n)2-S2-S2] 



(13) 
(14) 

(15) 



Given the Hamiltonian, the dynamics of the system is 
described by the equations 

dp 

lit 



dt 
dSg 

dt 



(Vs„i?) X S„ 



(16) 
(17) 
(18) 



These equations define a Poisson-system for y 
(p,x, Si,S2)^, i.e., 

dt 



B{y)VH, 



(19) 



with 



B{y) = 



Bi{y) 



B2{y) 



I 




(20) 



(21) 



(22) 



/ 


Vo 



Slz 

.-Sly 



S2Z 

~S2y 

A numerical scheme which preserves this special struc- 
ture can be expected to have a benevolent long-time be- 
haviour similar to the symplectic case, e.g., [21]. 

Furthermore, the post- Newtonian equations are in fact 
a perturbed Kepler problem, the corrections to the clas- 
sical motion scaling with ^ . In the units G = c = 1 this 
scaling is encoded in the higher orders of i or p^ in the 

post-Newtonian terms. As g > 1 and p^ < 1 in most 
circumstances, one has a Hamiltonian of the form 

H = H^ + 5H, SH^l, (23) 

where the 'larger' part can be solved analytically. Keep- 
ing in mind the just mentioned properties of the post- 
Newtonian equations we now present the already known 
integration methods. 



i/so(p, X, Si, S2) = 2 ^ . (12) The Poisson integrator suggested by [25] is designed 

9'^ to exactly preserve the structure (19). Starting with the 



POISSON INTEGRATOR FOR THE 
POST-NEWTONIAN EQUATIONS 
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Hamiltonian (23), the rclativistic contribution 5H is first 
split into an orbital part ^fpN,orb and a spin term -ffso.ss- 
The main idea is now to further decompose the spin-orbit 
and spin-spin parts as 



with 



and 



with 



Hso - Hgo + + Hqq, 



^so = ^(Soff ■e,)(L-e,), 



ss 



^ss 



^ss 



^ss -5—' 



Si - Si S2 • S2 



2q3 



2g3 



3 3(Si.n)(S2.n) 



'SS 



4 _ 3(Si.n)(Si.n) , 3(83 • n)(S2 • n) 



ss 



2q3 



(24) 

(25) 

(26) 

(27) 
(28) 
(29) 
(30) 



The major achievement of [25] was to find analytical so- 
lutions </5goi Vss ^'-'^ ^'-'^ '^^ each spin related part. 
This said, a structure preserving integrator $so,ss for 
the spinning terms is obtained by setting 



so,ss = ^so ° rso ° "Pso ° "Pss ° 'Pss ° 'Pss ° 'Pss- 

(31) 

Thus, if one solves the flow (/jn of the Newtonian part 
analytically and uses a symplectic scheme to calculate the 
orbital relativistic contributions $pN,orb, one can finally 
combine the three flows (p-^^, $pN.orb and <E>so,ss to obtain 
a structure preserving flow. In Section 6 we will discuss 
how to best arrange the individual flows. 



5. TRANSFORMATION TO CANONICAL 
FORM 

Instead of directly preserving the Poisson struc- 
ture (19) we can choose another way: The Darhoux-Lie 
theorem states that for every Poisson system (19) one can 
find a transformation 



*(y), 



(32) 



such that the system in the coordinates z is locally canon- 
ical. There are two properties of the post-Newtonian 
equations which enable us to find such a transforma- 
tion in this case. Firstly, the positions and momenta are 
already in canonical form. Therefore, a transformation 
(32) only has to focus on the spin coordinates. Secondly, 



by multiplying the equations of motions of the spins (18) 
with the respective spin Sq, we see that 



ld£SJ_dS, _ 
2 dt " dt ""^^ 



(33) 



i.e., the length of the individual spins is a first inte- 
gral. These two observations make it surprisingly easy 
to achieve the transformation to symplectic form. 

From the constancy of the spin-length we see that two 
spin variables are redundant. The post-Newtonian sys- 
tem can therefore be described by iV = 10 variables. 
Because of this, [23] proposed the use of cylindrical co- 
ordinates for the spins. Accordingly, we set 



fPa COs(^a)^ 
Sa = "T-aXa PaSin(^a) 



(34) 



where Xa relates the length of an object's spin to the 
square of its mass. The conservation of the spin- 
length allows for the elimination of one of the variables 
{Pa, 4>a,£.a)- Thus, we Can express pa in terms of as 



(35) 



whereby the spin and thus the Hamiltonian only depend 
on 0Q and ^q. 

In order to deduce the equations of motion for the 
two independent variables, we observe that the follow- 
ing equalities hold true: 



dH dH dSaa^ dH dSay 
- — + 



dSa 



dSa 



''ay "-'"fa 

dH _ dHdSax dH dSay dH dSaz 



(36) 



(37) 



d£,a dSax dS,a dSay dSaz 
-^T^ = -Pa Sin((/)a) = -Say, (38) 



dSay 



Pa COs{(f)a) =^ Sax, 



dcjja 

Saz = Xaml£,a- 



(39) 
(40) 



For the sake of shorter notation, we assume w.l.o.g. that 
XaiTT-a — 1 Until the end of this section. 
Due to relation (40), we have 



d^_dSa^_ dH_ dH_ 
dt - dt ~ dSax 95,/"^' 



(41) 



where the second equality is simply the equation of mo- 
tion for the z-componcnt of the spin. Substituting Sax 
and Say with the help of equations (38) and (39), and 
then applying (36) we get 



^ _dH_dSax _ dH dSay 

dt dSax d(j)a dSay d(t)a 



We now consider the time-derivatives of the x- and y- 
componcnts. Taking into account the equations of mo- 
tion for these components, the derivatives with regard to 
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time are 
dH 



dH 



Sa 



dH _ dS, 

Jay ~ 



dS. 



dt 



dH dSay 

dt 



dt d4>a 



dt ' 
(43) 



d^a dt dlpa 



dt 
(44) 



We can multiply the first equation with 



dSg 



and the 



second with and substract the two equations. This 
leads to 



dSax dSay 
d4>a d^a 



dH dSay „ 

'Jaz 



dSay dSax 
d(j>a d£,a 

dH dS„ 



dSay d^a 



dSa 



dt 

Jay 



dH 

dSay 



dSa 



d^a 



dia 
dH dSax 
dSax d£,a 



(45) 



Calculating the partial derivatives of the spin compo- 
nents with regard to the new variables on the left hand 
side and some of the partial derivatives on the right hand 
side, equation (45) becomes 



dt 



dH dSgy 

dSay d£,a 



dH 



<a 



dH dS„ 



dSaz dSax dia 



-ia. (46) 



Keeping in mind that 
relation (37), we arrive at 



9?a 



1 and then taking use of 



dt 



dH 

Wa' 



(47) 



All in all, the post-Newtonian equations for the ten in- 
dependent variables z = {p,£^a,'X-,4'a) read 



dz 

'dt 



d 
dt 



X 

01 
V02/ 



/Vp\ 

d^ 



H, 



(48) 



which is to say that the system in the new variables is 
symplectic. What is more, the transformation is defined 
globally as it is nothing other than expressing the spins 
with constant length via cylindrical coordinates. As a 
consequence, a structure preserving algorithm for the 
post-Newtonian equations can be obtained by carrying 
out the global transformation to canonical form and then 
applying a symplectic integrator. 



6. SCHEMES BASED ON SPLITTING 

6.1. On splitting methods 

It is well known, e.g., [26], that, given a Hamiltonian of 
the form (23), an integrator which is split in this natural 



way has a smaller local error than a comparable scheme. 
More precisely, suppose we were given some second order 
method. We could apply it with a given step size h to 
the whole system (23), thus constructing the flow $H.h,- 
But we could also apply the numerical scheme only to 
the 'small' part 5H and combine this symmetrically with 
the flow ip^ of the first term in (23). This would yield 
the second order integrators 



SHJi 



(49) 



and 



Now, if we compared the local errors, we would get 

yH,h-^H,h\\^Oih'') (51) 

for the numerical scheme applied to the whole system, 
but 



|<^ff,h-$spiit,h|| = 0(<5/i3 



(52) 
(53) 



for the splitting schemes. From this we observe that split- 
ting can reduce a scheme's local error. 

To see which of the two splitting methods is the better 
option, we first notice that for post-Newtonian equations, 
the relativistic parts are non-separable, i.e. the Hamilto- 
nian cannot be splitted in the form 



i?(p,x) = r(p) + F(x) 



(54) 



Unfortunately, when a system is non-separable, symplec- 
tic schemes have to be implicit, see e.g., [21]. As a con- 
sequence, a splitting integrator of the form (50) has to 
solve a system of implicit equations twice per time step 
whereas a scheme of the form (49) leads to only one im- 
plicit system per step. Thus, the splitting $spiit./i can 
be expected to be more efficient than $spiit,/f Numerical 
experiments by [24] have confirmed this so that we will 
only consider splittings of the form (49) in the following. 



6.2. On composition methods 

The drawback of a splitting scheme is that -no matter 
if we choose (49) or (50)- it is of second order even if the 
numerical scheme for the SH part is of (much) higher 
order. This can be overcome by clever composition: If 
we divide the step size h into smaller intervals h ~ aih + 
4- a^h -t- ... and set for some second order method 

<I'2ndJi 

<&comp,/i = 'I'2nd,Qi/i ° $2nd,a2/! ° ^2nd,a3h O (55) 

the thus obtained scheme $comp,h will be of higher order, 
provided that the satisfy specific conditions, see, e.g.. 
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[21], chapter II. If the underlying second order scheme 
^2nd,h is symplectic, $comp,/i, as a composition of many 
symplectic operations, will be so, too. 

Let us briefly state another useful fact about the im- 
plementation of composition schemes: If we choose the 
second order basic method as $2ndji = $spiitji, we have 

*I'comp,/i = •■■ ° ^2nd,aih ° ^2nd,ai+ih ° ■•• 



,.0.... (56) 



In the last step we could 'merge' terms thanks to the 
group property 



'Ph°fs= fh+s 



(57) 

which is valid for every exact flow, thus reducing the 
numerical effort. This would not be possible if we chose 
^2nd,h — 'i'spiit./i iustcad and, consequently, we found 
another advantage of splitting (49) over splitting (50). 

One of the most popular composition methods is the 
statc-of-thc-art Suzuki composition, [27], 

^4th,h — '&2nd.Q/i ° <I'2nd,Q/i ° $2nd,^J/i ° ^2nd,ah ° 2nd.Q/i, 

(58) 



with 



4^ 



4t 



■43 



(59) 
(60) 



This yields a 4th order method which is symmetric, i.e. 



4th, /i 



4th, -/i, 



(61) 



whenever the underlying scheme is. After all the back- 
ground information on splitting and composition meth- 
ods, we are now in the position to present structure pre- 
serving integration schemes for the post-Newtonian equa- 
tions which have been considered so far. 



• With the work of the previous two subsections and 
Section 5, we construct a Poisson integrator as fol- 
lows: We use the midpoint rule to calculate the 
flow $pN,orb corresponding to the orbital relativis- 
tic contributions. Then, we use the flow corre- 
sponding to the spin related parts as given in (31) 
and its adjoint 



*so.ss = *ss ° *ls ° ° *SS ° *SO ° *SO ° *SO: 

(( 



and symmetrically combine them with $pn., 
the form 



^Poisson 
6H,h 



so.ss, 



o <I) 



PN.orbJi 



o $ 



SO,SS,-i- 



(63) 



lb m 



(64) 



to obtain a numerical flow for all the relativistic 
parts 6H in (23). This numerical flow is symmetric 
and of second order as is any flow constructed in 
this way, sec, e.g., [21], chapter V. Therefore, we 
can combine it with the exact flow of the Newtonian 
part as in (49) to obtain the second order scheme 



^Poisson 
spht,/i 



^ ^Poisson 



°VN.i 



(65) 



This said, we can apply Suzuki's composition (58) 
which yields the 4th order symmetric Poisson inte- 
grator 



<j)Poisson (j)Poisson ^ 



Poisson ^jjPoisson ^j^Poisson ^Poisson 
SH,ah SH,ph SH,ah SH.ah ' 



(66) 



• In order to construct a symplectic scheme, we first 
apply the transformation to canonical form of Sec- 
tion 5. The Hamiltonian in the new variables z is 
still of the form (23). As a consequence, we can pro- 
ceed along the lines of the two subsections above. 
Therefore, we apply the implicit midpoint rule to 
the whole relativistic contribution 5H . This second 
order method can then be combined with the an- 
alytical solution of the Kepler problem, leading to 
the symplectic second order splitting scheme 



sympl 
split, /i 



SH,h 



(67) 



6.3. Splitting schemes for post-Newtonian 
equations 

We will present a Poisson integrator in accordance with 
[25] as well as a symplectic splitting scheme. In both 
cases we will use the implicit midpoint rule, already pro- 
posed in [17], which for any differential equation (1) has 
the form 



Yn+l 



- Yn+l 



(62) 



It is of second order and preserves symmetry and sym- 
plecticity, see, e.g., [21]. 



Again, we take use of Suzukis composition and ar- 
rive at the integrator 

$7™P' ^ $>5y"ipl ^ ^sympl ^ ^sympl ^ ^sympl ^ ^sympl 



SH,ah 



SH,ah 



SH,ah 



(68) 



which is symplectic and of order 4. 



The nice ideas behind them and their mathematical 
bounty notwithstanding, the just presented structure 
preserving algorithms based on splitting methods are not 
very efficient: Even in the case we use the group prop- 
erty (57) to 'merge' terms as illustrated in (56) whenever 
it is possible, the symplectic integrator $44")^' is still a 
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composition of 11 flows, five of which can only be com- 
puted via the solution of ten-dimensional imphcit sys- 
tems. Using the Poisson integrator (66) instead, we also 
have to calculate the midpoint rule five times. As it is 
only applied to the orbital motion, the implicit systems 
are reduced to 6 dimensions. But for this we have to pay 
heavily because, taking everything together, we have to 
calculate 67 flows during one time step, most of which 
are related to the spin contributions and require the cal- 
culations of numerous rotations, see [25]. All these facts, 
which will be confirmed in the numerical experiments sec- 
tion below, make us look for a more efficient alternative 
to solve the post-Newtonian equations of motion. This 
is where Gauss collocation methods come into play. 

7. GAUSS RUNGE-KUTTA METHODS 

Gauss-Runge-Kutta methods arc in fact collocation 
methods. Therefore, we give some background concern- 
ing these schemes. 



Here, li(t) denote the Lagrange-polynomials of degree s, 
hit)^l[^^. (76) 

Depending on which set of stages < ci < ... < Cg < 1 is 
chosen, different collocation methods can be constructed. 
By setting 

= ^(1 + q), (77) 

with Ci being the roots of the Legendre-polynomial of de- 
gree s, one obtains a Gauss collocation method. The order 
of this methods is 0(/i^*), cf. [28], which is optimal in the 
sense that there are no other s-stage one-step methods 
that achieve a similar high order without further numer- 
ical ruse. In addition. Gauss collocation methods are 
symplcctic and time-reversible, as is proven in [21]. Due 
to all these properties, Gauss-Runge-Kutta methods are 
quite natural candidates for the solution of non-separable 
Hamiltonian systems. 



7.1. On collocation polynomials 

Given an interval [to, to + h], stages < ci < ... < 
Cs < 1, and an initial- value problem of the form (1), the 
polynomial u{t) of degree s, satisfying 

u{tQ) - yo, (69) 

u{to + c,h) = f{to + Cih,u(to + Cih)), i = 1, s, 

(70) 

is called a collocation polynomial. In order to solve an 
initial-value problem by collocation, one has to find the 
polynomial u{t) which satisfies the collocation conditions 
(69), (70). This gives an approximate solution of the 
initial value problem after a time step h by setting 

y(to + /i)coi :=u(io + /i). (71) 

A detailed introduction to collocation methods can be 
found in [28]. 

It can now readily be shown, e.g., [28], that a colloca- 
tion method is equivalent to an implicit s-stage Runge- 
Kutta scheme 

s 

yn+i^yn + hY,hf{Y,), (72) 

s 

Y. = y„ + /i^a,,/(Y,), (73) 

with coefficients chosen as 

ay = / l3{t)dt, (74) 

bj = / k{t)dt. (75) 
Jo 



7.2. Gauss collocation for post-Newtonian 
equations 

In order to employ Gauss Runge-Kutta methods in 
post-Newtonian simulations, we just have to conduct the 
transformation (32) of Section 5 and then apply a Gauss 
collocation scheme to the whole system in the new co- 
ordinates z. Doing so, we will have to solve the system 
of implicit equations (73) for the inner stage values Y.^ 
during each time step. This system has s ■ 10 dimensions. 
Contrary to the splitting schemes, we have to solve the 
system only once when calculating the step z„ Zn-t-i- 
Besides, we can drastically reduce the effort for the so- 
lution of the implicit system if we take account of the 
following. 



7.3. Starting approximations 

An implicit system has to be solved iteratively. Of 
course, the number of iterations necessary to obtain the 
solution depends on the distance between the starting 
guesses Y° and the final values Y^. All the better then, 
if there were a fast method to obtain guesses that are very 
close to the final values. This is possible for the Gauss 
collocations' implicit systems: Given the inner stage val- 
ues of the last step z„_i — !■ z„, Y^^**' ^^""^ , we set 

Note that this requires no additional function evaluation 
as JiY^^'^^^ t,tcp^ j^^g j^^j calculated in the previous 
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step anyway. If the coefRcients /3y satisfy 



(79) 



one has, e.g., [21], ehapter VIII, 
\\Y,-Y^\\=0{h'), 



..,s. 



(80) 



The above sphtting schemes, in contrast, miss any sim- 
ilarly good starting approximations. Referring the inter- 
ested reader to [29] where we have hsted the coefficients 
Ci, bi, aij and Pij for s = 2,3,4,6, we now move on to 
the numerical tests. 



8. NUMERICAL EXPERIMENTS 

All simulations for this work were run on a Core 2 Duo 
E6600 machine with 2.4GHz and 4GB RAM. The codes 
for the simulations have been written in C++. 

In this section we test and compare the following algo- 
rithms: 

• Transformation to canonical form combined with 
Gauss Runge-Kutta methods for s = 2,3,4. The 
corresponding schemes are denoted by Gauss2, 
GaussS, and Gauss4, respectively. 



The symplectic splitting scheme 
referred to as Symp. 



which will be 



• The Poisson integrator (66), abbreviated by Poiss. 

• The classical 4th order explicit Runge-Kutta 
scheme given by the tableau 



(81) 


















1/2 


1/2 











1/2 





1/2 








1 
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1/6 


1/3 


1/3 


1/6. 



Hereafter, this method will by denoted by RK4. 
• The explicit Cash-Karp Runge-Kutta scheme 
















1/5 












3/l0 


3/40 


9/40 








3/5 


3/lU 


-9/10 


6/5 






1 


-11/54 


5/2 


-70/27 


35/27 




ys 


1631/55296 


175/512 


575/i3824 


44275/110592 


253/4096 




37/378 





250/621 


125/594 


512/1771, 



(82) 

as proposed by [30] , which is of order 5 and will be 
abbreviated by CK5. 



As the most reasonable measure for the efficiency, we 
compare the CPU calculation times. The algorithms' 
accuracy is tested with the help of the relative error in 
the Hamiltonian 



AH = 



H{yn) - H{yo) 



Hiyo) 

and the relative error along the trajectory 



(83) 



\ 



N 

E 

1=1 



(84) 



Here, superscript i denotes a vector's ith component. Un- 
less stated otherwise, the 'exact' solution yex(0 '^^iH be 
given by an s = 6-stage Gauss Runge-Kutta scheme with 
a step size ft, = 0.1 applied to the system in canonical 
coordinates. 

The simulations are aborted due to poor accuracy as 
soon as the error in the energy exceeds the tolerance 



AH > IQ-^. 



(85) 



At first glance, it seems arbitrary to subject the integra- 
tors to such an upper limit on the energy error. But we 
will show now that such a bound is indeed necessary. 



8.1. On the importance of energy conservation 

Let us assume there was no upper limit on the error 
in the energy and we applied RK4 to the orbital test 
case. For different step size h, this would yield the energy 
errors as given in Fig. 1. Let us now further assume we 
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0.0001 

1e-006 



<1 
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200000 300000 
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FIG. 1: For the classical RK4 scheme applied with different 
step sizes h to the purely orbital test case, the error in the 
energy AH is plotted against integration time t. 

wanted to plot Poincare sections for this two-dimensional 
problem in order to investigate it for chaotic behaviour. 
For different h, we would obtain the sections plotted in 
Fig. 2. For large /i, these resemble chaotic rather than 
the correct quasiperiodic motion. 



9 



0.01 r 
0.005 - 

£ 




h=5 




h=20 



h=40 
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30 35 40 45 50 55 60 30 35 40 45 50 55 60 65 
X X 



FIG. 2: Poincare sections at y = 0, py > for the purely 
orbital test case obtained with RK4 and four different step 
sizes h. 



Applying Gauss3 with the large step size h = 40 in- 
stead, the energy is conserved and consequently the sec- 
tions are calculated correctly, cf. Fig 3. We have thus 



35 40 45 50 55"" 




The other relevant parameter is the factor Xa, already 
introduced in section 5, that links masses with spins via 



|Sa|| = Xaml. 



(89) 



Hence, the nature of a binary's orbit depends on the pa- 
rameters a,xi,X2 and the initial values 

z(0) = 

(p,(0),py(0),p,(0), ei(0), 6(0), xiO), y(0), z(0), 0i(O), 02(0))' 

(90) 

This said, the three kinds of motion are represented by 
the following respective examples: 

• With the set of initial data 



z(0) 



0,— ,0,0,0,35,0,0,0,0 



(91) 



Xi = X2 = 0, 



we model a system without spin effects. The 
spin contributions being switched off, the post- 
Newtonian system is integrable, e.g., [23], and the 
motion is restricted to the initial plane due to 
the conservation of the angular momentum. We 
present the orbit and the Poincare sections for 
t € [0, lO"^] as obtained via 'exact' integration in 
Fig. 4. The motion is apparently quasipcriodic. 



FIG. 3: The left panel shows Poincare sections at y = 0, 
Py > for the purely orbital test case obtained with GaussS 
and h = 40. In the right panel, the corresponding error in the 
energy AH is plotted against integration time t. 

illustrated that a threshold for the relative error in the 
energy is inevitable if wc want to obtain reliable infor- 
mation on the chaoticity. Let us now present the test 
cases with the help of which we compare the individual 
methods. 



40 45 So 55 60 B5 



FIG. 4: For the test case without spin contributions and t G 
[0, 10^], the left panel shows an extract of the trajectory. The 
Poincare sections for y = and py > are given in the right 
panel. 



8.2. The test cases 

We model three different kinds of motion, each of which 
is often encountered in binary simulations. Wc always fix 
the total mass as m = 1. Consequently, the important 
parameter concerning the two compact object's masses 
is the mass ratio cr = — . The individual masses and the 
reduced mass are thus given as 



• As a second test case, we choose the data set 

y(0) = (^0, ^,0, 0.25, -0.025, 35, 0, 0, J, , 
1 

3' 



Xi = X2 



(92) 



mi 



m2 



A* = 



" 1 + ct' 
1 

a 



(86) 
(87) 
(88) 



In Fig. 5, we plot a part of the orbital trajectory for 
t G [0,10'']. Alongside this, we plot the frequency 
spectrum of the x component for Ii = [0, 10^] and 
h = [10^ - 10^, 10^]. We see that although the spin 
contributions have been switched on, the motion is 
still regular. 



10 





































































tonn 


m 




114 (1. 


elm 







FIG. 5: For the test case (92), the left panels shows the tra- 
jectory for t £ [0,50000]. The frequency spectra \f^{iu)\ for 
the time intervals h = [0, 10^] and h = [lO'^ - lO", lO'^] are 
depicted in the right panel. 



• We also consider a chaotic orbit. More precisely, 
we set 



y(0) = ( 1,0,^,0,0.25,-0.025,6,0,0,^,^ 



cr = 1, 



(93) 



We illustrate the chaotic behaviour by showing a 
part of the orbital trajectory and the FLI in Fig. 6. 
The FLI shows characteristically chaotic traits, cf. 
[31]. 
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FIG. 6: For the chaotic test case, the left panel shows the 
trajectory for t £ [0,25000]. The linearly growing FLI is 
depicted in the right panel in semi-logarithmic scale. 

Having thus established the test cases, we are able to 
start with our experiments. 



8.3. Comparing the splitting schemes 

We first compare the two splitting schemes. As they 
are exactly the same in the non-spinning case, we turn 
towards the regular spinning example (92) and plot the 
respective error in the Hamiltonian for various h in Fig. 7. 
We see no difference in the accuracy. But when compar- 
ing the corresponding calculation times in table I we see 
that the Poisson scheme is much slower. 

Testing the splitting schemes for the chaotic test case 
we see that Symp falls victim to criterion 85 for step sizes 
as small as /i = 5 but it can cope with it for /i < 1. Not 
so Poiss which even fails for h = 0.01. Consequently, 
the symplectic splitting (68) is superior to the Poisson 
integrator (66). But as we will corroborate now, it is by 
now means the best option for post-Newtonian systems. 




^1e-013 
<l 

1e-014 



4e+006 8B+006 2e+006 6e+006 

t t 



FIG. 7: For initial data (92), t G [0,10''] and different step 
sizes h, the relative error in the Hamiltonian AH is plotted 
against time t for the splitting integrators of section ??. No 
difference can be spotted between them. 



Integrator 


/!. = 40/i = 20 h = 5 h = l 


Symp 
Poiss 


9.90 18.72 67.75 304.95 
17.37 34.23 133.34 655.11 



TABLE I: The CPU calculation times in [s] for the two split- 
ting integrators applied to the regular, spinning test case (92) 
with different step sizes h. The integration interval was 
t G [0, 10^] in all simulations. 



8.4. Comparing integration schemes 

Here, we compare the symplectic splitting to the (ex- 
plicit and structure preserving) Runge-Kutta schemes. 
First, we list the calculation times for simulations with 
the orbital test case in table II. As would have been 
expected, the explicit schemes are faster than the other 
methods for equal step sizes. But they have to be applied 
with small step sizes in order not to hurt the constraint 
on the energy error. We also see that Symp is by far the 
slowest algorithm. Doing the same observations for the 



Integrator 


h = AQ 


/i = 20 


= 5 


h = 1 


h = 0.5 


h = 0.1 


RK4 


a 


a 


a 


13.80 


27.58 


137.91 


CK5 


a 


a 


4.70 


23.01 


46.07 


230.01 


Gauss2 


a 


3.89 


11.44 


43.81 


81.73 


344.41 


Gauss3 


3.27 


5.32 


15.44 


58.48 


105.73 


422.69 


Gauss4 


3.96 


6.47 


18.74 


67.26 


120.59 


443.49 


Symp 


4.36 


8.35 


30.95 


142.86 







TABLE II: The CPU calculation times in [s] for several 
schemes applied to the orbital test case (92) with different 
step sizes h. The integration interval was t G [0, 10''] in all 
simulations, 'a' signifies 'aborted due to condition (85)'. 

regular spinning case, we get equal results, cf. table HI. 

As the errors of the individual integrators behave sim- 
ilarly for both regular orbits, we only show the case 
with spins included. We plot the error along the trajec- 
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Integrator 


h = 40 


ft = 20 h = 5 


h = 1 


h = 0.5 


h = 0.1 


RK4 


a 


a 


a 


43.99 


87.93 


439.44 


CK5 


a 


a 


14.31 


71.56 


143.02 


716.16 


Gauss2 


a 


10.03 


29.37 


111.23 


205.68 


852.88 


GaussS 


8.10 


13.19 


38.09 


141.19 


255.60 


997.65 


Gauss4 


10.00 


16.44 


46.63 


160.95 


283.26 


1068.56 


Symp 


9.90 


18.72 


67.75 


304.95 







TABLE III: The CPU calculation times in [s] for several 
schemes applied to the regular, spinning test case (92) with 
different step sizes h. The integration interval was t € [0, 10^] 
in all simulations, 'a' signifies 'aborted due to condition (85)'. 



arising from the initial conditions (93). Again, we start 
with hsting the calculation times of simulations with var- 
ious step sizes in table (IV) , right after which we plot the 
error along the trajectory in Fig. 10 and the relative error 
in the energy for various simulations in Fig. 11. The first 
point to mention here is that the explicit methods require 
prohibitively small step sizes in order not to exceed the 
error bar (85). As of the structure preserving candidates, 
the result is qualitatively the same as in the regular sim- 
ulations: Symp seems to be better than Gauss2 which 
struggles with the chaotic case. But it obviously cannot 
match the performance of the fast and accurate Gauss3 
and Gauss4. 



tory (84) in Fig. 8 and the relative error in the Hamilto- 
nian (83) in Fig. 9. We see that although Symp is more 
accurate than CK5 and Gauss2 with equal step sizes, 
it has a larger error than the much faster GaussS and 
Gauss4 with equal or even much larger step sizes. 
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FIG. 8: For initial data (92) and t G [0,10'^], the relative 
error along the trajectory, cf. (84), is plotted against time t 
for various integration schemes. 
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FIG. 9: For initial data (92) and t G [0, lO'^], the relative error 
in the Hamiltonian Aff is plotted against time t for various 
integration schemes. 

We now turn our attention towards the chaotic motion 



Integrator 


/i = 5 


h= 1 


h = 0.5 


ft = 0.1 ft 


= 0.05 ft = 0.01 


RK4 


a 


a 


a 


a 


a 2997.76 


CK5 


a 


a 


a 


a 


a 4840.96 


Gauss2 


a 


a 


a 


1190.14 




Gauss3 


a 


a 


449.47 


1548.56 




Gauss4 


a 


347.22 


566.67 


1893.45 




Symp 


a 


463.78 


833.85 


3445.49 





TABLE IV: The CPU calculation times in [s] for several 
schemes applied to the chaotic test case (92) with different 
step sizes ft. The integration interval was t G [0, 10^] in all 
simulations, 'a' signifies 'aborted due to condition (85)'. 
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FIG. 10: For initial data (93), the relative error along the 
trajectory, cf. (84), is plotted against time t for various inte- 
gration schemes. 

One interesting point which stands out for all three 
initial data is that the difference in CPU times between 
explicit and Gauss Runge-Kutta scheme decreases for 
smaller step sizes. This is thanks to the starting approx- 
imations introduced in subsection 7.3. The smaller the 
step size, the closer the initial guess of the iterations gets 
to the correct values due to relation (80). Consequently 
the average number of iterations per step decreases along- 
side h. To illustrate this, we list the iterations per step 
of Gauss4 in table V. 

We have seen that the structure preserving algorithms 
have excellent conservation properties when applied to 
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FIG. 11: For initial data (93) and t £ [0,10''], the relative 
error in the Hamiltonian AH is plotted against time t for 
various integration schemes. 



test case ft = 40 /i = 20 ft = 5 /i = 1 /i = 0.5 /i = 0.1 

initial values (91) 9.19 fM 5l6 3M 2^99 2l3~ 

initial values (92) 9.31 7.59 5.30 3.54 3.07 2.24 

initial values (93) a 9.21 7.44 4.85 



TABLE V: Number of iterations per step for Gauss4, applied 
with different step sizes to the three test cases, 'a' signifies 
'aborted due to condition (85)'. 
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FIG. 12: The radial distance q as function of integration time 
t for the regular spinning orbit with radiation efi'ects included. 



In the subsections above, CK5 and Gauss3/Gauss4 
showed the best results for explicit and structure pre- 
serving schemes, respectively. We thus focus on these 
integrators and compare their performance with the ra- 
diation turned on. We first list the calculation times for 
the three schemes applied with different step sizes each, 
cf. table VI. With increasing time steps, the differ- 
ence in CPU time becomes ever smaller as the collocation 
methods' average number of iterations per step decreases 
analogously to the conservative case. 

As a measure for the accuracy we plot the relative error 
along the trajectory (84) in Fig. 13. Taking into account 
the calculation times, the collocation methods yield the 
better results for less computational costs - just as in the 
conservative case. 



symplectic systems. What will happen if we add a radi- 
ation term to the binary system? 

8.5. Systems with radiation 

Adding a dissipative term, the system loses the struc- 
ture which gave rise to the advantageous integrators in 
the first place. But it is known from classical mechan- 
ics that, at least in this field, structure preserving al- 
gorithms outperform explicit schemes also when a non- 
conservative term is added to the Hamiltonian. In order 
to examine the corresponding behaviour for relativistic 
binaries, we restrict ourselves to the initial data (92) and 
modify the equation of motion of the momenta (16) to 
account for radiation. We choose a model for the radia- 
tion force Frad derived by [32] which is commonly used 
in general relativity and set 

5 = -Vxi? + F,.ad. (94) 

at 

To illustrate its effects on the trajectory, we plot the evo- 
lution of the radial distance q for our regular, spinning 
test case (92) as given by the exact solution in Fig. 12. 
Here, we calculate the 'exact' solution with CK5 and the 
very small step size /i = 0.01. As time increases, the 
distance between the two particles is decreasing faster 
and faster. For t > 500 000, the post-Newtonian approx- 
imation will soon lose its validity. Thus, we restrict our 
simulations to an interval t G [0, 500 000]. 
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FIG. 13: The relative error along the trajectory, err against 
integration time t for explicit and collocation schemes. 



Integrator 


ft = 40 ft = 20 ft = 5 ft = 1 ft = 0.5 ft = 0.1 


CK5 
Gauss3 
Gauss4 


0.14 0.28 1.21 5.75 11.50 57.51 
0.92 1.47 4.04 16.45 25.69 99.77 
1.08 1.72 4.72 17.28 32.00 133.36 



TABLE VI: The CPU calculation times in [s] for explicit and 
implicit Runge-Kutta schemes applied with different step sizes 
ft to the test case (92) with radiation effects included. The 
integration interval was t £ [0, 500 000] for all simulations. 
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9. SUMMARY 

We have seen that structure preserving algorithms are 
necessary for the long-time integration of post-Newtonian 
equations of motion as they guarantee the conservation of 
the energy which is inevitable in investigations for chaos. 
Thus, in this work we analysed several algorithms - a 
Poisson integrator based on the Poisson structure, a sym- 
plectic splitting scheme and Gauss Runge-Kutta meth- 
ods. We observed large discrepancies in the performance 
of the individual structure preserving methods. Some 
even fared worse than explicit methods. More specifi- 
cally, the Poisson integrator turned out to be extremely 
slow when applied to our test cases. The symplectic 
scheme based on state-of-the-art splitting and compo- 
sition techniques could compete with a Gauss Runge- 
Kutta scheme with two inner stages but was completely 
out-beaten by Gauss collocation schemes with three or 



more inner stages. These Gauss methods turned out to 
be by far the most efficient and most accurate option. 
Even for dissipative systems, they delivered more accu- 
rate results for equal computational cost than high order 
explicit Rungc-Kutta schemes. Therefore, we strongly 
recommend to use a transformation of the system to sym- 
plectic form combined with a Gauss Runge-Kutta scheme 
for the numerical long-time analysis of post-Newtonian 
systems. 
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